Data and code for paper on evolution during a 2017 epidemic of Pasteuria ramosa in Little Appleton Lake. Authors: Camden D. Gowler Haley Essington Patrick A. Clay Bruce O’Brien Clara L. Shaw Rebecca Bilich Meghan A. Duffy

Analysis of data collected on Little Appleton’s 2017 Pasteuria epidemic

Field data on infection prevalence and host density

lilApp <- readr::read_csv(here("data/lilApp.csv"))

## New names:
## * `` -> ...1

## Rows: 11 Columns: 54

## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (3): Date, Lake, Counted.By
## dbl (51): ...1, Total, X, UA, UJ, Uninfected.Males, Uninfected.Ephip, A.Past...

## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# do this to see which parasites are there (past, micg, spiro)
ttt <- lilApp %>% 
  select(Uninfected.Males:M.Past...Spider) %>%
  gather("Parasite", "Count",Uninfected.Males:M.Past...Spider) %>%
  group_by(Parasite) %>%
  summarise(totinfe = sum(Count)) %>%
  filter(totinfe > 0)

past <- lilApp %>% # find columns with pasteuria in title
  select(contains("past"))
lilApp$pasteuria.inf <- rowSums(past)

micg <- lilApp %>% 
  select(contains("mic"))
lilApp$micg.inf <- rowSums(micg)

spiro <- lilApp %>% 
  select(contains("spiro"))
lilApp$spiro.inf <- rowSums(spiro)

LilApp_allpara <- lilApp %>%
  mutate(pasteuria.prev = pasteuria.inf/Total,
         micg.prev = micg.inf/Total,
         spiro.prev = spiro.inf/Total)

# super contrived way of doing this 
sample_vec <- c(rep(NA, 11))
collect <- c(3, 5, 7)
for (i in 1:length(collect)) {
  x <- collect[i]
  sample_vec[x] <- paste0("collected",i)
}
sample_vec[is.na(sample_vec)] <- "regular"
LilApp_allpara$sampling <- sample_vec
LilApp_allpara$sampling <- as.factor(LilApp_allpara$sampling)

str(LilApp_allpara)

## spec_tbl_df [11 × 61] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ ...1                     : num [1:11] 1 2 3 4 5 6 7 8 9 10 ...
##  $ Date                     : chr [1:11] "7/27/17" "8/10/17" "8/28/17" "9/8/17" ...
##  $ Lake                     : chr [1:11] "LittleAppleton" "LittleAppleton" "LittleAppleton" "LittleAppleton" ...
##  $ Counted.By               : chr [1:11] "MAR" "KKH" "MAR" "KKH" ...
##  $ Total                    : num [1:11] 265 308 359 321 206 249 116 171 203 154 ...
##  $ X                        : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ UA                       : num [1:11] 70 97 49 46 22 101 36 62 90 41 ...
##  $ UJ                       : num [1:11] 185 179 183 163 71 51 30 83 100 88 ...
##  $ Uninfected.Males         : num [1:11] 0 0 0 0 0 2 1 13 12 22 ...
##  $ Uninfected.Ephip         : num [1:11] 0 0 0 0 0 0 0 0 0 3 ...
##  $ A.Pasteuria              : num [1:11] 0 7 6 17 18 3 5 3 1 0 ...
##  $ J.Pasteuria              : num [1:11] 0 10 61 82 55 79 31 5 0 0 ...
##  $ Male.Pasteuria           : num [1:11] 0 0 0 0 0 1 0 0 0 0 ...
##  $ Ephip.Pastueria          : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ Brood                    : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ E.Brood                  : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ A.MicG                   : num [1:11] 8 13 39 12 19 11 10 5 0 0 ...
##  $ J.MicG                   : num [1:11] 0 0 20 1 13 1 0 0 0 0 ...
##  $ Male.MicG                : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ Ephip.MicG               : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ A.Spiro                  : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ J.Spiro                  : num [1:11] 2 2 1 0 0 0 0 0 0 0 ...
##  $ Male.Spiro               : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ Ephip.Spiro              : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ A.White.Spiro            : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ J.White.Spiro            : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ Male.White.Spiro         : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ E.White.Spiro            : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ A.Yellow.Spiro           : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ A.Metsch                 : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ J.Metsch                 : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ Male.Metsch              : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ Ephip.Metsch             : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ A.Larssonia              : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ A.Spiderman              : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ J.Spiderman              : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ E.Spiderman              : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ Male.Spiderman           : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ A.MicG...Past            : num [1:11] 0 0 0 0 4 0 1 0 0 0 ...
##  $ J.MicG...Past            : num [1:11] 0 0 0 0 4 0 2 0 0 0 ...
##  $ A.MicG...Spiro           : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ A.Metsch...Brood         : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ AMicG...Brood            : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ A.MicG...Metsch          : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ J.MicG...Metsch          : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ M.MicG...Metsch          : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ A.Metsch...Past          : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ J.Metsch...Past          : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ A.MicG...White.Spiro     : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ Male.white.spiro...metsch: num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ A.Past...Spider          : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ M.Past...Spider          : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
##  $ Julian.Day               : num [1:11] 208 222 240 251 261 276 280 286 300 304 ...
##  $ Density                  : num [1:11] 235212 305732 107071 128571 67081 ...
##  $ pasteuria.inf            : num [1:11] 0 17 67 99 81 83 39 8 1 0 ...
##  $ micg.inf                 : num [1:11] 8 13 59 13 40 12 13 5 0 0 ...
##  $ spiro.inf                : num [1:11] 2 2 1 0 0 0 0 0 0 0 ...
##  $ pasteuria.prev           : num [1:11] 0 0.0552 0.1866 0.3084 0.3932 ...
##  $ micg.prev                : num [1:11] 0.0302 0.0422 0.1643 0.0405 0.1942 ...
##  $ spiro.prev               : num [1:11] 0.00755 0.00649 0.00279 0 0 ...
##  $ sampling                 : Factor w/ 4 levels "collected1","collected2",..: 4 4 1 4 2 4 3 4 4 4 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   ...1 = col_double(),
##   ..   Date = col_character(),
##   ..   Lake = col_character(),
##   ..   Counted.By = col_character(),
##   ..   Total = col_double(),
##   ..   X = col_double(),
##   ..   UA = col_double(),
##   ..   UJ = col_double(),
##   ..   Uninfected.Males = col_double(),
##   ..   Uninfected.Ephip = col_double(),
##   ..   A.Pasteuria = col_double(),
##   ..   J.Pasteuria = col_double(),
##   ..   Male.Pasteuria = col_double(),
##   ..   Ephip.Pastueria = col_double(),
##   ..   Brood = col_double(),
##   ..   E.Brood = col_double(),
##   ..   A.MicG = col_double(),
##   ..   J.MicG = col_double(),
##   ..   Male.MicG = col_double(),
##   ..   Ephip.MicG = col_double(),
##   ..   A.Spiro = col_double(),
##   ..   J.Spiro = col_double(),
##   ..   Male.Spiro = col_double(),
##   ..   Ephip.Spiro = col_double(),
##   ..   A.White.Spiro = col_double(),
##   ..   J.White.Spiro = col_double(),
##   ..   Male.White.Spiro = col_double(),
##   ..   E.White.Spiro = col_double(),
##   ..   A.Yellow.Spiro = col_double(),
##   ..   A.Metsch = col_double(),
##   ..   J.Metsch = col_double(),
##   ..   Male.Metsch = col_double(),
##   ..   Ephip.Metsch = col_double(),
##   ..   A.Larssonia = col_double(),
##   ..   A.Spiderman = col_double(),
##   ..   J.Spiderman = col_double(),
##   ..   E.Spiderman = col_double(),
##   ..   Male.Spiderman = col_double(),
##   ..   A.MicG...Past = col_double(),
##   ..   J.MicG...Past = col_double(),
##   ..   A.MicG...Spiro = col_double(),
##   ..   A.Metsch...Brood = col_double(),
##   ..   AMicG...Brood = col_double(),
##   ..   A.MicG...Metsch = col_double(),
##   ..   J.MicG...Metsch = col_double(),
##   ..   M.MicG...Metsch = col_double(),
##   ..   A.Metsch...Past = col_double(),
##   ..   J.Metsch...Past = col_double(),
##   ..   A.MicG...White.Spiro = col_double(),
##   ..   Male.white.spiro...metsch = col_double(),
##   ..   A.Past...Spider = col_double(),
##   ..   M.Past...Spider = col_double(),
##   ..   Julian.Day = col_double(),
##   ..   Density = col_double()
##   .. )
##  - attr(*, "problems")=<externalptr>

summary(LilApp_allpara)

##       ...1          Date               Lake            Counted.By       
##  Min.   : 1.0   Length:11          Length:11          Length:11         
##  1st Qu.: 3.5   Class :character   Class :character   Class :character  
##  Median : 6.0   Mode  :character   Mode  :character   Mode  :character  
##  Mean   : 6.0                                                           
##  3rd Qu.: 8.5                                                           
##  Max.   :11.0                                                           
##      Total             X           UA               UJ        Uninfected.Males
##  Min.   :116.0   Min.   :0   Min.   : 22.00   Min.   : 30.0   Min.   : 0.000  
##  1st Qu.:187.0   1st Qu.:0   1st Qu.: 43.50   1st Qu.: 75.5   1st Qu.: 0.000  
##  Median :220.0   Median :0   Median : 62.00   Median : 88.0   Median : 1.000  
##  Mean   :233.8   Mean   :0   Mean   : 62.73   Mean   :110.3   Mean   : 6.636  
##  3rd Qu.:286.5   3rd Qu.:0   3rd Qu.: 83.00   3rd Qu.:171.0   3rd Qu.:12.500  
##  Max.   :359.0   Max.   :0   Max.   :101.00   Max.   :185.0   Max.   :23.000  
##  Uninfected.Ephip  A.Pasteuria      J.Pasteuria    Male.Pasteuria   
##  Min.   : 0.000   Min.   : 0.000   Min.   : 0.00   Min.   :0.00000  
##  1st Qu.: 0.000   1st Qu.: 0.500   1st Qu.: 0.00   1st Qu.:0.00000  
##  Median : 0.000   Median : 3.000   Median :10.00   Median :0.00000  
##  Mean   : 3.909   Mean   : 5.455   Mean   :29.36   Mean   :0.09091  
##  3rd Qu.: 0.000   3rd Qu.: 6.500   3rd Qu.:58.00   3rd Qu.:0.00000  
##  Max.   :40.000   Max.   :18.000   Max.   :82.00   Max.   :1.00000  
##  Ephip.Pastueria     Brood      E.Brood      A.MicG          J.MicG      
##  Min.   :0       Min.   :0   Min.   :0   Min.   : 0.00   Min.   : 0.000  
##  1st Qu.:0       1st Qu.:0   1st Qu.:0   1st Qu.: 3.00   1st Qu.: 0.000  
##  Median :0       Median :0   Median :0   Median :10.00   Median : 0.000  
##  Mean   :0       Mean   :0   Mean   :0   Mean   :10.73   Mean   : 3.182  
##  3rd Qu.:0       3rd Qu.:0   3rd Qu.:0   3rd Qu.:12.50   3rd Qu.: 1.000  
##  Max.   :0       Max.   :0   Max.   :0   Max.   :39.00   Max.   :20.000  
##    Male.MicG   Ephip.MicG    A.Spiro     J.Spiro         Male.Spiro
##  Min.   :0   Min.   :0    Min.   :0   Min.   :0.0000   Min.   :0   
##  1st Qu.:0   1st Qu.:0    1st Qu.:0   1st Qu.:0.0000   1st Qu.:0   
##  Median :0   Median :0    Median :0   Median :0.0000   Median :0   
##  Mean   :0   Mean   :0    Mean   :0   Mean   :0.4545   Mean   :0   
##  3rd Qu.:0   3rd Qu.:0    3rd Qu.:0   3rd Qu.:0.5000   3rd Qu.:0   
##  Max.   :0   Max.   :0    Max.   :0   Max.   :2.0000   Max.   :0   
##   Ephip.Spiro A.White.Spiro J.White.Spiro Male.White.Spiro E.White.Spiro
##  Min.   :0    Min.   :0     Min.   :0     Min.   :0        Min.   :0    
##  1st Qu.:0    1st Qu.:0     1st Qu.:0     1st Qu.:0        1st Qu.:0    
##  Median :0    Median :0     Median :0     Median :0        Median :0    
##  Mean   :0    Mean   :0     Mean   :0     Mean   :0        Mean   :0    
##  3rd Qu.:0    3rd Qu.:0     3rd Qu.:0     3rd Qu.:0        3rd Qu.:0    
##  Max.   :0    Max.   :0     Max.   :0     Max.   :0        Max.   :0    
##  A.Yellow.Spiro    A.Metsch    J.Metsch  Male.Metsch  Ephip.Metsch  A.Larssonia
##  Min.   :0      Min.   :0   Min.   :0   Min.   :0    Min.   :0     Min.   :0   
##  1st Qu.:0      1st Qu.:0   1st Qu.:0   1st Qu.:0    1st Qu.:0     1st Qu.:0   
##  Median :0      Median :0   Median :0   Median :0    Median :0     Median :0   
##  Mean   :0      Mean   :0   Mean   :0   Mean   :0    Mean   :0     Mean   :0   
##  3rd Qu.:0      3rd Qu.:0   3rd Qu.:0   3rd Qu.:0    3rd Qu.:0     3rd Qu.:0   
##  Max.   :0      Max.   :0   Max.   :0   Max.   :0    Max.   :0     Max.   :0   
##   A.Spiderman  J.Spiderman  E.Spiderman Male.Spiderman A.MicG...Past   
##  Min.   :0    Min.   :0    Min.   :0    Min.   :0      Min.   :0.0000  
##  1st Qu.:0    1st Qu.:0    1st Qu.:0    1st Qu.:0      1st Qu.:0.0000  
##  Median :0    Median :0    Median :0    Median :0      Median :0.0000  
##  Mean   :0    Mean   :0    Mean   :0    Mean   :0      Mean   :0.4545  
##  3rd Qu.:0    3rd Qu.:0    3rd Qu.:0    3rd Qu.:0      3rd Qu.:0.0000  
##  Max.   :0    Max.   :0    Max.   :0    Max.   :0      Max.   :4.0000  
##  J.MicG...Past    A.MicG...Spiro A.Metsch...Brood AMicG...Brood A.MicG...Metsch
##  Min.   :0.0000   Min.   :0      Min.   :0        Min.   :0     Min.   :0      
##  1st Qu.:0.0000   1st Qu.:0      1st Qu.:0        1st Qu.:0     1st Qu.:0      
##  Median :0.0000   Median :0      Median :0        Median :0     Median :0      
##  Mean   :0.5455   Mean   :0      Mean   :0        Mean   :0     Mean   :0      
##  3rd Qu.:0.0000   3rd Qu.:0      3rd Qu.:0        3rd Qu.:0     3rd Qu.:0      
##  Max.   :4.0000   Max.   :0      Max.   :0        Max.   :0     Max.   :0      
##  J.MicG...Metsch M.MicG...Metsch A.Metsch...Past J.Metsch...Past
##  Min.   :0       Min.   :0       Min.   :0       Min.   :0      
##  1st Qu.:0       1st Qu.:0       1st Qu.:0       1st Qu.:0      
##  Median :0       Median :0       Median :0       Median :0      
##  Mean   :0       Mean   :0       Mean   :0       Mean   :0      
##  3rd Qu.:0       3rd Qu.:0       3rd Qu.:0       3rd Qu.:0      
##  Max.   :0       Max.   :0       Max.   :0       Max.   :0      
##  A.MicG...White.Spiro Male.white.spiro...metsch A.Past...Spider M.Past...Spider
##  Min.   :0            Min.   :0                 Min.   :0       Min.   :0      
##  1st Qu.:0            1st Qu.:0                 1st Qu.:0       1st Qu.:0      
##  Median :0            Median :0                 Median :0       Median :0      
##  Mean   :0            Mean   :0                 Mean   :0       Mean   :0      
##  3rd Qu.:0            3rd Qu.:0                 3rd Qu.:0       3rd Qu.:0      
##  Max.   :0            Max.   :0                 Max.   :0       Max.   :0      
##    Julian.Day       Density       pasteuria.inf      micg.inf    
##  Min.   :208.0   Min.   :  4013   Min.   : 0.00   Min.   : 0.00  
##  1st Qu.:245.5   1st Qu.:  8457   1st Qu.: 0.50   1st Qu.: 3.00  
##  Median :276.0   Median : 14333   Median :17.00   Median :12.00  
##  Mean   :267.7   Mean   : 81362   Mean   :35.91   Mean   :14.91  
##  3rd Qu.:293.0   3rd Qu.:117821   3rd Qu.:74.00   3rd Qu.:13.00  
##  Max.   :317.0   Max.   :305732   Max.   :99.00   Max.   :59.00  
##    spiro.inf      pasteuria.prev       micg.prev         spiro.prev      
##  Min.   :0.0000   Min.   :0.000000   Min.   :0.00000   Min.   :0.000000  
##  1st Qu.:0.0000   1st Qu.:0.002463   1st Qu.:0.01689   1st Qu.:0.000000  
##  Median :0.0000   Median :0.055195   Median :0.04050   Median :0.000000  
##  Mean   :0.4545   Mean   :0.151335   Mean   :0.06050   Mean   :0.001530  
##  3rd Qu.:0.5000   3rd Qu.:0.320872   3rd Qu.:0.08013   3rd Qu.:0.001393  
##  Max.   :2.0000   Max.   :0.393204   Max.   :0.19417   Max.   :0.007547  
##        sampling
##  collected1:1  
##  collected2:1  
##  collected3:1  
##  regular   :8  
##                
##

LA_pastprev_plot <- ggplot() +
  geom_line(data = LilApp_allpara, aes(x = Julian.Day, y = pasteuria.prev*100), size = 1.1, alpha = 0.7) +
  geom_point(data = LilApp_allpara, aes(x = Julian.Day, y = pasteuria.prev*100, fill=sampling), size = 5, alpha = 1.0, shape = 21, show.legend = FALSE) +
  ylim(0, 0.45*100) +
  scale_x_continuous(breaks = c(213, 244, 274, 305), labels = c("Aug-1", "Sept-1","Oct-1", "Nov-1")) +
  scale_fill_manual(values=c("#5445b1", "#749dae", "#f3c483",'#7E7E7E')) +
  xlab("") +
  ylab("Prevalence \n(%)") 
LA_pastprev_plot

LA_density_plot <- ggplot() +
  geom_line(data = LilApp_allpara, aes(x = Julian.Day, y = Density), size = 1.1, alpha = 0.7) +
  geom_point(data = LilApp_allpara, aes(x = Julian.Day, y = Density, fill=sampling), size = 5, alpha = 1.0, shape = 21, show.legend = FALSE) +
  scale_y_log10() +
  scale_x_continuous(breaks = c(213, 244, 274, 305), labels = c("Aug-1", "Sept-1","Oct-1", "Nov-1")) +
  scale_fill_manual(values=c("#5445b1", "#749dae", "#f3c483",'#7E7E7E')) +
  xlab("Date") +
  ylab("*D. dentifera* density \n(per m^2)") +
  theme(axis.title.y = ggtext::element_markdown())
LA_density_plot

fielddataplot <- plot_grid(LA_pastprev_plot, LA_density_plot, labels = "auto", ncol = 1, align = "v")
fielddataplot

ggsave(here("figures", "fielddataplot.jpg"), fielddataplot, units = "in", width = 5, height = 7, dpi = 300)

Infection Assay Data

Checking out sample sizes

# Load the data for analyzing virulence
LA_infected <- readr::read_csv(here("data/LA_infected.csv"))

## Rows: 267 Columns: 13

## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (4): rep, treatment, infection, para_time_point
## dbl  (7): clone, block, host_time_point, lifespan, avg, spores_tot, clutches
## date (2): death_date, birth_date

## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# Get rid of time 4 (which had only one host clone)
LA_infected <- subset(LA_infected, host_time_point != "4")

# Calculating the proportion that never reproduced
LA_infected <- LA_infected %>%
  mutate(anyrepro = if_else(clutches == "0", "0", "1"))

LA_infected$anyrepro <- as.numeric(LA_infected$anyrepro)

# Adding a column for log spores
LA_infected$logspores <- log10(LA_infected$spores_tot)
LA_infected$lnspores <- log(LA_infected$spores_tot)

LA_infected_samplesizes <- LA_infected %>%
  count(clone, para_time_point)

# save
write.csv(LA_infected_samplesizes, here("tables", "LA_infected_samplesizes.csv"), row.names = F)

Infected vs. uninfected comparison

LA_infected$clone <- as.factor(as.character(LA_infected$clone))
LA_infected$host_time_point <- as.factor(as.character(LA_infected$host_time_point))
LA_infected$para_time_point <- as.factor(as.character(LA_infected$para_time_point))

LA_infected$exposed <- ifelse(LA_infected$para_time_point == "none", 'unexposed', 'infected')

# Note: all of the animals in this dataset are either controls or infected, so don't need to worry about exposed but uninfected animals

LA_infected_summary <- LA_infected %>%
  group_by(exposed, host_time_point, clone, para_time_point) %>%
  summarise(meanlifespan = mean(lifespan), meanclutches = mean(clutches))

## `summarise()` has grouped output by 'exposed', 'host_time_point', 'clone'. You can override using the `.groups` argument.

lifespaninfvcontrol <- ggplot(LA_infected_summary, 
                               aes(x=exposed, y=meanlifespan)) + 
  geom_violin() + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, color = '#5c1a33', size = 2) + 
  ylab("Mean lifespan \n(days)")  +
  theme_cowplot() +
  theme(axis.title.x = element_blank())


lifespaninfvcontrol

# making a plot comparing exposed vs. unexposed
fecundityinfvcontrol <- ggplot(LA_infected_summary, 
                               aes(x=exposed, y=meanclutches)) + 
  geom_violin() + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, color = '#5c1a33', size = 2) +
  ylab("Mean # clutches")  +
  theme_cowplot() +
  theme(axis.title.x = element_blank())

fecundityinfvcontrol

Let’s arrange the lifespan & fecundity plots into a single figure:

infvcontrolplot <- plot_grid(lifespaninfvcontrol, fecundityinfvcontrol, labels = "auto", ncol = 1, align = "v")
infvcontrolplot

ggsave(here("figures", "infvcontrolplot.jpg"), infvcontrolplot, units = "in", width = 5, height = 8, dpi = 300)

infvcontrollifespanmodel<-glmer(lifespan ~ exposed + (1 | clone), family = "poisson", data=LA_infected) 

summary(infvcontrollifespanmodel)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: lifespan ~ exposed + (1 | clone)
##    Data: LA_infected
## 
##      AIC      BIC   logLik deviance df.resid 
##   2047.5   2058.2  -1020.8   2041.5      256 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.7830 -0.9877  0.1476  0.9803  4.7174 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  clone  (Intercept) 0.003289 0.05735 
## Number of obs: 259, groups:  clone, 28
## 
## Fixed effects:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       3.53313    0.01679  210.41   <2e-16 ***
## exposedunexposed  0.36282    0.02649   13.69   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## exposdnxpsd -0.353

# Testing for overdispersion (using code from Ben Bolker via Michelle Fearon)
overdisp_fun <- function(model) {
  rdf <- df.residual(model)
  rp <- residuals(model,type="pearson")
  Pearson.chisq <- sum(rp^2)
  prat <- Pearson.chisq/rdf
  pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
  c(chisq=Pearson.chisq,ratio=prat,rdf=rdf,p=pval)
}

overdisp_fun(infvcontrollifespanmodel)

##        chisq        ratio          rdf            p 
## 5.605749e+02 2.189746e+00 2.560000e+02 7.929284e-25

#very overdispersed, so let's go to a negative binomial distribution
infvcontrollifespanmodel.nb <-glmer.nb(lifespan ~ exposed + (1 | clone), data=LA_infected) 

overdisp_fun(infvcontrollifespanmodel.nb)

##       chisq       ratio         rdf           p 
## 231.3587907   0.9072894 255.0000000   0.8534423

summary(infvcontrollifespanmodel.nb)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Negative Binomial(23.7957)  ( log )
## Formula: lifespan ~ exposed + (1 | clone)
##    Data: LA_infected
## 
##      AIC      BIC   logLik deviance df.resid 
##   1913.6   1927.9   -952.8   1905.6      255 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3581 -0.5694  0.0871  0.6945  2.8266 
## 
## Random effects:
##  Groups Name        Variance  Std.Dev.
##  clone  (Intercept) 7.208e-05 0.00849 
## Number of obs: 259, groups:  clone, 28
## 
## Fixed effects:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       3.53253    0.01835 192.540   <2e-16 ***
## exposedunexposed  0.36737    0.04151   8.849   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## exposdnxpsd -0.439

infvcontrolrepro<-glmer(clutches ~ exposed + (1 | clone), family = "poisson", data=LA_infected) 

overdisp_fun(infvcontrolrepro)

##        chisq        ratio          rdf            p 
## 8.038865e+02 3.140182e+00 2.560000e+02 7.119246e-58

infvcontrolrepro.nb <-glmer.nb(clutches ~ exposed + (1 | clone), data=LA_infected) 

overdisp_fun(infvcontrolrepro.nb)

##        chisq        ratio          rdf            p 
## 303.84415988   1.19154573 255.00000000   0.01934465

summary(infvcontrolrepro.nb)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Negative Binomial(0.8023)  ( log )
## Formula: clutches ~ exposed + (1 | clone)
##    Data: LA_infected
## 
##      AIC      BIC   logLik deviance df.resid 
##    836.9    851.1   -414.5    828.9      255 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.8719 -0.5923 -0.4879  0.2711  6.3781 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  clone  (Intercept) 0.466    0.6827  
## Number of obs: 259, groups:  clone, 28
## 
## Fixed effects:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -0.4457     0.1873  -2.379   0.0174 *  
## exposedunexposed   2.8101     0.2343  11.995   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## exposdnxpsd -0.281

Infection prevalence in infection assays

# The dataset we were working with above only has infected animals or controls. Let's upload the dataset that will also include animals that were exposed but did not become infected
LA_full <- readr::read_csv(here("data/LA_full.csv"))

## Rows: 374 Columns: 13

## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (4): rep, treatment, infection, para_time_point
## dbl  (7): clone, block, host_time_point, lifespan, avg, spores_tot, clutches
## date (2): death_date, birth_date

## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# since we want to look at infectivity, let's get rid of the control/unexposed animals
LA_exposed <- LA_full %>%
  subset(para_time_point != "none")

# get rid of time 4 hosts, since we didn't have parasites from that time point
LA_exposed <- LA_exposed %>%
  subset(host_time_point != '4')

# filter out ones that are missing spore counts
LA_exposed <- LA_exposed %>%
  drop_na(spores_tot)

# Looking for clones that had very few animals in the experiment
LA_exposed_samplesizes <- LA_exposed %>%
  count(clone, para_time_point)

LA_exposed_samplesizes$remove <- ifelse(LA_exposed_samplesizes$n < 3, 'y', 'n')

# save
write.csv(LA_exposed_samplesizes, here("tables", "LA_exposed_samplesizes.csv"), row.names = F)


# Using spore yield to determine if it was infected; this is more sensitive than the 'infection' column, which was based on a quick visual assessment at death -- that inspection missed some that had very low spore densities
LA_exposed$infected <- ifelse(LA_exposed$spores_tot == 0, 0, 1)


# Calculate number infected and total sample sizes
LA_exposed_prevalence <- LA_exposed %>%
  group_by(clone, host_time_point, para_time_point) %>%
  summarise(total = n(), inf = sum(infected))

## `summarise()` has grouped output by 'clone', 'host_time_point'. You can override using the `.groups` argument.

# Calculate proportion infected per clone
LA_exposed_prevalence$propinf <- LA_exposed_prevalence$inf/LA_exposed_prevalence$total

LA_exposed_prevalence$host_time_point <- as.character(LA_exposed_prevalence$host_time_point)

infprevplot <- ggplot(LA_exposed_prevalence, 
                     aes(x=host_time_point, y=propinf, fill=para_time_point)) + 
  geom_violin(position=position_dodge(1)) + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, show.legend = FALSE) +
  scale_fill_manual(values=c("#5445b1", "#749dae", "#f3c483")) + 
  scale_color_manual(values=c("#5445b1", "#749dae", "#f3c483")) +
  xlab("Host time point") +ylab("Proportion infected") +
  labs(color= "Parasite time point") + labs(fill="Parasite time point") +
  theme_cowplot()

infprevplot

infprevplota <- LA_exposed_prevalence %>%
  filter(host_time_point == para_time_point) %>%
  ggplot(aes(x=host_time_point, y=propinf, fill=para_time_point)) + 
  geom_violin(position=position_dodge(1), show.legend = FALSE) + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, show.legend = FALSE) +
  scale_fill_manual(values=c("#5445b1", "#749dae", "#f3c483")) + 
  scale_color_manual(values=c("#5445b1", "#749dae", "#f3c483")) +
  ylim(0,1.00) +
  xlab("Host time point") +ylab("Proportion infected (contemporary parasites)") +
  labs(color= "Parasite time point") + labs(fill="Parasite time point") +
  theme_cowplot()

infprevplota

infprevplotb <- LA_exposed_prevalence %>%
  filter(host_time_point == 3) %>%
  ggplot(aes(x=para_time_point, y=propinf, fill=para_time_point)) + 
  geom_violin(position=position_dodge(1), show.legend = FALSE) + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, show.legend = FALSE) +
  scale_fill_manual(values=c("#C0C0C0", "#f3c483")) + 
  scale_color_manual(values=c("#C0C0C0", "#f3c483")) +
  ylim(0,1.00) +
  xlab("Parasite time point") +ylab("Proportion infected (time 3 hosts)") +
  theme_cowplot()

infprevplotb

infprevplot2 <- plot_grid(infprevplota, infprevplotb, labels = "auto", ncol = 2, align = "v")
infprevplot2

ggsave(here("figures", "infprevplot.jpg"), infprevplot2, units = "in", width = 8, height = 5, dpi = 300)

LA_exposed_prevalence_sametime <- subset(LA_exposed_prevalence, host_time_point == para_time_point)

prevmodel1<-glmer(cbind(inf, total-inf) ~ para_time_point + (1 | clone), family='binomial', data=LA_exposed_prevalence_sametime) 

# Testing for overdispersion (using code from Ben Bolker via Michelle Fearon)
overdisp_fun <- function(model) {
  rdf <- df.residual(model)
  rp <- residuals(model,type="pearson")
  Pearson.chisq <- sum(rp^2)
  prat <- Pearson.chisq/rdf
  pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
  c(chisq=Pearson.chisq,ratio=prat,rdf=rdf,p=pval)
}

overdisp_fun(prevmodel1)

##      chisq      ratio        rdf          p 
## 15.3492797  0.6673600 23.0000000  0.8816798

summary(prevmodel1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(inf, total - inf) ~ para_time_point + (1 | clone)
##    Data: LA_exposed_prevalence_sametime
## 
##      AIC      BIC   logLik deviance df.resid 
##    102.2    107.4    -47.1     94.2       23 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.60835 -0.54861 -0.08292  0.75747  1.13735 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  clone  (Intercept) 0.3203   0.5659  
## Number of obs: 27, groups:  clone, 27
## 
## Fixed effects:
##                  Estimate Std. Error z value Pr(>|z|)
## (Intercept)        0.1359     0.4972   0.273    0.785
## para_time_point2   0.4872     0.5938   0.820    0.412
## para_time_point3   0.8869     0.5818   1.524    0.127
## 
## Correlation of Fixed Effects:
##             (Intr) pr_t_2
## par_tm_pnt2 -0.851       
## par_tm_pnt3 -0.873  0.750

LA_exposed_prevalence_sametime$para_time_point <- as.factor(LA_exposed_prevalence_sametime$para_time_point)

# Testing for overall effect of parasite time point
drop1(prevmodel1,test="Chisq")

## Single term deletions
## 
## Model:
## cbind(inf, total - inf) ~ para_time_point + (1 | clone)
##                 npar    AIC    LRT Pr(Chi)
## <none>               102.17               
## para_time_point    2 100.88 2.7111  0.2578

# Now comparing the two sets of spores that were used to infect time 3 hosts
LA_exposed_prevalence_time3 <- subset(LA_exposed_prevalence, host_time_point == '3')

prevmodel2<-glmer(cbind(inf, total-inf) ~ para_time_point + (1 | clone), family='binomial', data=LA_exposed_prevalence_time3) 

overdisp_fun(prevmodel2)

##      chisq      ratio        rdf          p 
## 27.7920278  1.1580012 24.0000000  0.2689096

summary(prevmodel2)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(inf, total - inf) ~ para_time_point + (1 | clone)
##    Data: LA_exposed_prevalence_time3
## 
##      AIC      BIC   logLik deviance df.resid 
##     94.4     98.3    -44.2     88.4       24 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3300 -0.8751  0.2081  0.8352  1.4208 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  clone  (Intercept) 0.1248   0.3532  
## Number of obs: 27, groups:  clone, 14
## 
## Fixed effects:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        1.1495     0.2654   4.331 1.49e-05 ***
## para_time_point3  -0.1808     0.3297  -0.548    0.584    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## par_tm_pnt3 -0.673

# Getting chi-square to be parallel to initial analysis
drop1(prevmodel2,test="Chisq")

## Single term deletions
## 
## Model:
## cbind(inf, total - inf) ~ para_time_point + (1 | clone)
##                 npar    AIC     LRT Pr(Chi)
## <none>               94.402                
## para_time_point    1 92.703 0.30031  0.5837

Contemporary virulence

Comparing hosts exposed to parasites from the same time point

LA_sametime <- LA_infected %>%
  filter(para_time_point != "none") 

LA_sametime$para_time_point <- as.factor(LA_sametime$para_time_point)
LA_sametime$host_time_point <- as.factor(as.character(LA_sametime$host_time_point))

LA_sametime$para_time_point <- droplevels(LA_sametime$para_time_point)

LA_sametime <- subset(LA_sametime, host_time_point == para_time_point)

LA_sametime$paragrowth <- (LA_sametime$spores_tot)/(LA_sametime$lifespan)


LA_sametime_summary <- LA_sametime %>%
  group_by(clone, para_time_point) %>%
  summarise(meanlifespan = mean(lifespan), meanclutches = mean(clutches), meanspores = mean(spores_tot), meanlogspores = mean(logspores), meanparagrowth = mean(paragrowth))

## `summarise()` has grouped output by 'clone'. You can override using the `.groups` argument.

contemplifespanplot <- ggplot(LA_sametime_summary, 
                     aes(x=para_time_point, y=meanlifespan, fill=para_time_point)) + 
  geom_violin(show.legend = FALSE) + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0, seed = 1), alpha = 0.5, show.legend = FALSE) +
  scale_fill_manual(values=c("#5445b1", "#749dae", "#f3c483")) + 
  scale_color_manual(values=c("#5445b1", "#749dae", "#f3c483")) +
  xlab("Parasite time point") +ylab("Mean lifespan \n(days)") +
  labs(color= "Parasite time point") + labs(fill="Parasite time point") +
  theme_cowplot()


contemplifespanplot

contemporarylifespan <-glmer(lifespan ~ host_time_point + (1 | clone), family = "poisson", data=LA_sametime) 

overdisp_fun(contemporarylifespan)

##        chisq        ratio          rdf            p 
## 2.543648e+02 1.829963e+00 1.390000e+02 8.579296e-09

summary(contemporarylifespan)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: lifespan ~ host_time_point + (1 | clone)
##    Data: LA_sametime
## 
##      AIC      BIC   logLik deviance df.resid 
##   1080.7   1092.5   -536.3   1072.7      139 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0520 -0.9930 -0.0005  0.9847  4.0506 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  clone  (Intercept) 0.009985 0.09992 
## Number of obs: 143, groups:  clone, 27
## 
## Fixed effects:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       3.49051    0.07526  46.376   <2e-16 ***
## host_time_point2  0.09306    0.08537   1.090    0.276    
## host_time_point3  0.06512    0.08299   0.785    0.433    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) hst__2
## hst_tm_pnt2 -0.878       
## hst_tm_pnt3 -0.910  0.799

contemporarylifespan.nb <-glmer.nb(lifespan ~ host_time_point + (1 | clone), data=LA_sametime) 

overdisp_fun(contemporarylifespan.nb)

##       chisq       ratio         rdf           p 
## 136.4245942   0.9885840 138.0000000   0.5219351

summary(contemporarylifespan.nb)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Negative Binomial(28.2922)  ( log )
## Formula: lifespan ~ host_time_point + (1 | clone)
##    Data: LA_sametime
## 
##      AIC      BIC   logLik deviance df.resid 
##   1038.2   1053.0   -514.1   1028.2      138 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9149 -0.7669  0.0226  0.6781  3.6005 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  clone  (Intercept) 0.001694 0.04116 
## Number of obs: 143, groups:  clone, 27
## 
## Fixed effects:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       3.41497    0.07576  45.077   <2e-16 ***
## host_time_point2  0.15898    0.08331   1.908   0.0563 .  
## host_time_point3  0.14369    0.08314   1.728   0.0840 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) hst__2
## hst_tm_pnt2 -0.885       
## hst_tm_pnt3 -0.924  0.817

drop1(contemporarylifespan.nb,test="Chisq")

## Single term deletions
## 
## Model:
## lifespan ~ host_time_point + (1 | clone)
##                 npar    AIC    LRT Pr(Chi)
## <none>               1038.2               
## host_time_point    2 1037.1 2.9334  0.2307

emmeans(contemporarylifespan.nb, specs = pairwise ~ host_time_point)

## $emmeans
##  host_time_point emmean     SE  df asymp.LCL asymp.UCL
##  1                 3.41 0.0758 Inf      3.27      3.56
##  2                 3.57 0.0389 Inf      3.50      3.65
##  3                 3.56 0.0317 Inf      3.50      3.62
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast estimate     SE  df z.ratio p.value
##  1 - 2     -0.1590 0.0833 Inf  -1.908  0.1363
##  1 - 3     -0.1437 0.0831 Inf  -1.728  0.1947
##  2 - 3      0.0153 0.0504 Inf   0.304  0.9504
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates

contempreproplot <- ggplot(LA_sametime_summary, 
                     aes(x=para_time_point, y=meanclutches, fill=para_time_point)) + 
  geom_violin(show.legend = FALSE) + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0, seed = 1), alpha = 0.5, show.legend = FALSE) +
  scale_fill_manual(values=c("#5445b1", "#749dae", "#f3c483")) + 
  scale_color_manual(values=c("#5445b1", "#749dae", "#f3c483")) +
  xlab("Parasite time point") +ylab("Mean number of clutches") +
  labs(color= "Parasite time point") + labs(fill="Parasite time point") +
  theme_cowplot()


contempreproplot

contemporaryclutches <-glmer(clutches ~ host_time_point + (1 | clone), family = "poisson", data=LA_sametime) 

overdisp_fun(contemporaryclutches)

##        chisq        ratio          rdf            p 
## 2.724757e+02 1.960257e+00 1.390000e+02 1.038315e-10

summary(contemporaryclutches)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: clutches ~ host_time_point + (1 | clone)
##    Data: LA_sametime
## 
##      AIC      BIC   logLik deviance df.resid 
##    445.8    457.6   -218.9    437.8      139 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2287 -0.7444 -0.4875  0.0935  4.9567 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  clone  (Intercept) 1.627    1.275   
## Number of obs: 143, groups:  clone, 27
## 
## Fixed effects:
##                  Estimate Std. Error z value Pr(>|z|)
## (Intercept)       -0.1837     0.7839  -0.234    0.815
## host_time_point2  -0.1733     0.8959  -0.193    0.847
## host_time_point3  -0.7804     0.8713  -0.896    0.370
## 
## Correlation of Fixed Effects:
##             (Intr) hst__2
## hst_tm_pnt2 -0.858       
## hst_tm_pnt3 -0.875  0.760

contemporaryclutches.nb <-glmer.nb(clutches ~ host_time_point + (1 | clone), data=LA_sametime) 

overdisp_fun(contemporaryclutches.nb)

##       chisq       ratio         rdf           p 
## 106.5188352   0.7718756 138.0000000   0.9783568

summary(contemporaryclutches.nb)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Negative Binomial(0.3527)  ( log )
## Formula: clutches ~ host_time_point + (1 | clone)
##    Data: LA_sametime
## 
##      AIC      BIC   logLik deviance df.resid 
##    344.8    359.6   -167.4    334.8      138 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.5582 -0.4579 -0.4028  0.1347  3.7312 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  clone  (Intercept) 0.8264   0.9091  
## Number of obs: 143, groups:  clone, 27
## 
## Fixed effects:
##                  Estimate Std. Error z value Pr(>|z|)
## (Intercept)        0.2162     0.7690   0.281    0.779
## host_time_point2  -0.4843     0.8720  -0.555    0.579
## host_time_point3  -1.0601     0.8472  -1.251    0.211
## 
## Correlation of Fixed Effects:
##             (Intr) hst__2
## hst_tm_pnt2 -0.863       
## hst_tm_pnt3 -0.884  0.784

drop1(contemporaryclutches.nb,test="Chisq")

## Single term deletions
## 
## Model:
## clutches ~ host_time_point + (1 | clone)
##                 npar    AIC    LRT Pr(Chi)
## <none>               344.82               
## host_time_point    2 342.78 1.9565   0.376

emmeans(contemporaryclutches.nb, specs = pairwise ~ host_time_point)

## $emmeans
##  host_time_point emmean    SE  df asymp.LCL asymp.UCL
##  1                0.216 0.769 Inf     -1.29    1.7234
##  2               -0.268 0.441 Inf     -1.13    0.5965
##  3               -0.844 0.397 Inf     -1.62   -0.0658
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast estimate    SE  df z.ratio p.value
##  1 - 2       0.484 0.872 Inf   0.555  0.8437
##  1 - 3       1.060 0.847 Inf   1.251  0.4229
##  2 - 3       0.576 0.566 Inf   1.017  0.5658
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates

contempsporesplot <- ggplot(LA_sametime_summary, 
                     aes(x=para_time_point, y=meanspores, fill=para_time_point)) + 
  geom_violin(show.legend = FALSE) + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0, seed = 1), alpha = 0.5, show.legend = FALSE) +
  scale_fill_manual(values=c("#5445b1", "#749dae", "#f3c483")) + 
  scale_color_manual(values=c("#5445b1", "#749dae", "#f3c483")) +
  xlab("Parasite time point") +ylab("Mean number of \nspores per infected host") +
  labs(color= "Parasite time point") + labs(fill="Parasite time point") +
  scale_y_log10() +
  theme_cowplot()


contempsporesplot

contemporarylogspores <- lmer(logspores ~ host_time_point + (1 | clone), data=LA_sametime) 

overdisp_fun(contemporarylogspores)

##       chisq       ratio         rdf           p 
##  48.7586242   0.3533234 138.0000000   1.0000000

summary(contemporarylogspores)

## Linear mixed model fit by REML ['lmerMod']
## Formula: logspores ~ host_time_point + (1 | clone)
##    Data: LA_sametime
## 
## REML criterion at convergence: 284.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8804 -0.3749  0.1072  0.6637  1.5019 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  clone    (Intercept) 0.05924  0.2434  
##  Residual             0.37508  0.6124  
## Number of obs: 143, groups:  clone, 27
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)       4.96189    0.21368  23.221
## host_time_point2 -0.07965    0.24515  -0.325
## host_time_point3 -0.25471    0.23558  -1.081
## 
## Correlation of Fixed Effects:
##             (Intr) hst__2
## hst_tm_pnt2 -0.872       
## hst_tm_pnt3 -0.907  0.791

drop1(contemporarylogspores,test="Chisq")

## Single term deletions
## 
## Model:
## logspores ~ host_time_point + (1 | clone)
##                 npar    AIC   LRT Pr(Chi)
## <none>               287.95              
## host_time_point    2 286.11 2.166  0.3386

emmeans(contemporarylogspores, specs = pairwise ~ host_time_point)

## $emmeans
##  host_time_point emmean     SE   df lower.CL upper.CL
##  1                 4.96 0.2201 18.9     4.50     5.42
##  2                 4.88 0.1209 22.2     4.63     5.13
##  3                 4.71 0.0995 20.1     4.50     4.91
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast estimate    SE   df t.ratio p.value
##  1 - 2      0.0796 0.251 19.6   0.317  0.9462
##  1 - 3      0.2547 0.242 19.1   1.055  0.5527
##  2 - 3      0.1751 0.157 21.3   1.118  0.5136
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates

#R3 asked about the non-log-transformed data, so let's do that
contemporaryspores <- lmer(spores_tot ~ host_time_point + (1 | clone), data=LA_sametime) 

overdisp_fun(contemporaryspores)

##        chisq        ratio          rdf            p 
## 3.097726e+12 2.244729e+10 1.380000e+02 0.000000e+00

summary(contemporaryspores)

## Linear mixed model fit by REML ['lmerMod']
## Formula: spores_tot ~ host_time_point + (1 | clone)
##    Data: LA_sametime
## 
## REML criterion at convergence: 3760.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.4286 -0.5585 -0.3070  0.3365  4.7327 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  clone    (Intercept) 2.312e+09  48086  
##  Residual             2.339e+10 152932  
## Number of obs: 143, groups:  clone, 27
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)        136891      47944   2.855
## host_time_point2    35071      55183   0.636
## host_time_point3   -22837      52941  -0.431
## 
## Correlation of Fixed Effects:
##             (Intr) hst__2
## hst_tm_pnt2 -0.869       
## hst_tm_pnt3 -0.906  0.787

drop1(contemporaryspores,test="Chisq")

## Single term deletions
## 
## Model:
## spores_tot ~ host_time_point + (1 | clone)
##                 npar    AIC   LRT Pr(Chi)
## <none>               3837.7              
## host_time_point    2 3836.6 2.914  0.2329

contempparagrowthplot <- ggplot(LA_sametime_summary, 
                     aes(x=para_time_point, y=meanparagrowth, fill=para_time_point)) + 
  geom_violin(show.legend = FALSE) + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0, seed = 1), alpha = 0.5, show.legend = FALSE) +
  scale_fill_manual(values=c("#5445b1", "#749dae", "#f3c483")) + 
  scale_color_manual(values=c("#5445b1", "#749dae", "#f3c483")) +
  xlab("Parasite time point") +ylab("Parasite growth rate \n(spores per day)") +
  labs(color= "Parasite time point") + labs(fill="Parasite time point") +
  theme_cowplot()


contempparagrowthplot

contemporaryparasitegrowth <- lmer(paragrowth ~ host_time_point + (1 | clone), data=LA_sametime) 

overdisp_fun(contemporaryparasitegrowth)

##      chisq      ratio        rdf          p 
## 1994481979   14452768        138          0

summary(contemporaryparasitegrowth)

## Linear mixed model fit by REML ['lmerMod']
## Formula: paragrowth ~ host_time_point + (1 | clone)
##    Data: LA_sametime
## 
## REML criterion at convergence: 2736
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6309 -0.5728 -0.2677  0.3091  4.2680 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  clone    (Intercept)  2065795 1437    
##  Residual             15245387 3905    
## Number of obs: 143, groups:  clone, 27
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)        4279.9     1313.2   3.259
## host_time_point2    538.4     1507.8   0.357
## host_time_point3  -1219.3     1448.3  -0.842
## 
## Correlation of Fixed Effects:
##             (Intr) hst__2
## hst_tm_pnt2 -0.871       
## hst_tm_pnt3 -0.907  0.790

drop1(contemporaryparasitegrowth,test="Chisq")

## Single term deletions
## 
## Model:
## paragrowth ~ host_time_point + (1 | clone)
##                 npar    AIC    LRT Pr(Chi)
## <none>               2791.7               
## host_time_point    2 2791.4 3.6498  0.1612

emmeans(contemporaryparasitegrowth, specs = pairwise ~ host_time_point)

## $emmeans
##  host_time_point emmean   SE   df lower.CL upper.CL
##  1                 4280 1357 17.4     1422     7138
##  2                 4818  746 22.0     3272     6365
##  3                 3061  613 19.9     1782     4339
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast estimate   SE   df t.ratio p.value
##  1 - 2        -538 1548 18.4  -0.348  0.9358
##  1 - 3        1219 1489 17.8   0.819  0.6965
##  2 - 3        1758  965 21.1   1.821  0.1870
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates

Parasite evolution comparison

Comparing hosts from time 3 that were exposed to parasites from time 1 vs. time 3

LA_paraevol <- LA_infected %>%
  filter(host_time_point == "3") %>%
  filter(para_time_point != "none")

LA_paraevol$paragrowth <- (LA_paraevol$spores_tot)/(LA_paraevol$lifespan)


LA_paraevol_summary <- LA_paraevol %>%
  group_by(clone, para_time_point) %>%
  summarise(meanlifespan = mean(lifespan), meanclutches = mean(clutches), meanspores = mean(spores_tot), meanlogspores = mean(logspores), meanparagrowth = mean(paragrowth))

## `summarise()` has grouped output by 'clone'. You can override using the `.groups` argument.

Lifespan

lifespanparaevol <- ggplot(LA_paraevol_summary, 
                     aes(x=para_time_point, y=meanlifespan, fill=para_time_point)) + 
  geom_violin(show.legend = FALSE) + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0, seed = 1), alpha = 0.5, show.legend = FALSE) +
  scale_fill_manual(values=c("#C0C0C0", "#f3c483")) + 
  scale_color_manual(values=c("#C0C0C0", "#f3c483")) +
  xlab("Parasite time point") +ylab("Mean lifespan \n(days)") +
  labs(color= "Parasite time point") + labs(fill="Parasite time point") +
  theme_cowplot()


lifespanparaevol

paraevollifespan <-glmer(lifespan ~ para_time_point + (1 | clone), family = "poisson", data=LA_paraevol) 

overdisp_fun(paraevollifespan)

##        chisq        ratio          rdf            p 
## 2.459747e+02 1.732216e+00 1.420000e+02 1.428517e-07

summary(paraevollifespan)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: lifespan ~ para_time_point + (1 | clone)
##    Data: LA_paraevol
## 
##      AIC      BIC   logLik deviance df.resid 
##   1057.3   1066.2   -525.7   1051.3      142 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.2121 -0.8940  0.1315  0.8178  3.1930 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  clone  (Intercept) 0.003669 0.06057 
## Number of obs: 145, groups:  clone, 14
## 
## Fixed effects:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       3.50151    0.02748 127.421   <2e-16 ***
## para_time_point3  0.05260    0.02914   1.805   0.0711 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## par_tm_pnt3 -0.577

paraevollifespan.nb <-glmer.nb(lifespan ~ para_time_point + (1 | clone), data=LA_paraevol) 

overdisp_fun(paraevollifespan.nb)

##       chisq       ratio         rdf           p 
## 137.3167343   0.9738775 141.0000000   0.5720026

summary(paraevollifespan.nb)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Negative Binomial(38.342)  ( log )
## Formula: lifespan ~ para_time_point + (1 | clone)
##    Data: LA_paraevol
## 
##      AIC      BIC   logLik deviance df.resid 
##   1027.4   1039.3   -509.7   1019.4      141 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1564 -0.6213  0.0274  0.6221  2.4806 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  clone  (Intercept) 0.001195 0.03458 
## Number of obs: 145, groups:  clone, 14
## 
## Fixed effects:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       3.50845    0.03044 115.264   <2e-16 ***
## para_time_point3  0.04913    0.03937   1.248    0.212    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## par_tm_pnt3 -0.680

Reproduction

reproparaevol <- ggplot(LA_paraevol_summary, 
                     aes(x=para_time_point, y=meanclutches, fill=para_time_point)) + 
  geom_violin(show.legend = FALSE) + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0, seed = 1), alpha = 0.5, show.legend = FALSE) +
  scale_fill_manual(values=c("#C0C0C0", "#f3c483")) + 
  scale_color_manual(values=c("#C0C0C0", "#f3c483")) +
  xlab("Parasite time point") +ylab("Mean number of clutches") +
  labs(color= "Parasite time point") + labs(fill="Parasite time point") +
  theme_cowplot()


reproparaevol

paraevolclutches <-glmer(clutches ~ para_time_point + (1 | clone), family = "poisson", data=LA_paraevol) 

overdisp_fun(paraevolclutches)

##        chisq        ratio          rdf            p 
## 3.081137e+02 2.169815e+00 1.420000e+02 2.583562e-14

summary(paraevolclutches)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: clutches ~ para_time_point + (1 | clone)
##    Data: LA_paraevol
## 
##      AIC      BIC   logLik deviance df.resid 
##    344.7    353.6   -169.3    338.7      142 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5607 -0.7308 -0.5406 -0.2792  7.3191 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  clone  (Intercept) 0.5547   0.7448  
## Number of obs: 145, groups:  clone, 14
## 
## Fixed effects:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -1.0469     0.2896  -3.615 0.000301 ***
## para_time_point3   0.4432     0.2380   1.862 0.062623 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## par_tm_pnt3 -0.516

paraevolclutches.nb <-glmer.nb(clutches ~ para_time_point + (1 | clone), data=LA_paraevol)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00817007 (tol = 0.002, component 1)

## boundary (singular) fit: see ?isSingular

overdisp_fun(paraevolclutches.nb)

##       chisq       ratio         rdf           p 
## 158.6689373   1.1253116 141.0000000   0.1467461

summary(paraevolclutches.nb)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Negative Binomial(0.2552)  ( log )
## Formula: clutches ~ para_time_point + (1 | clone)
##    Data: LA_paraevol
## 
##      AIC      BIC   logLik deviance df.resid 
##    276.5    288.4   -134.3    268.5      141 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.4292 -0.4292 -0.3989  0.2190  7.3489 
## 
## Random effects:
##  Groups Name        Variance  Std.Dev. 
##  clone  (Intercept) 2.695e-10 1.642e-05
## Number of obs: 145, groups:  clone, 14
## 
## Fixed effects:
##                  Estimate Std. Error z value Pr(>|z|)   
## (Intercept)       -0.8615     0.2975  -2.895  0.00379 **
## para_time_point3   0.4492     0.4023   1.117  0.26419   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## par_tm_pnt3 -0.739
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular

Spores

sporesparaevol <- ggplot(LA_paraevol_summary, 
                     aes(x=para_time_point, y=meanspores, fill=para_time_point)) + 
  geom_violin(show.legend = FALSE) + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0, seed = 1), alpha = 0.5, show.legend = FALSE) +
  scale_fill_manual(values=c("#C0C0C0", "#f3c483")) + 
  scale_color_manual(values=c("#C0C0C0", "#f3c483")) +
  xlab("Parasite time point") +ylab("Mean number of \nspores per infected host") +
  labs(color= "Parasite time point") + labs(fill="Parasite time point") +
  scale_y_log10() +
  theme_cowplot()


sporesparaevol

paraevollogspores <- lmer(logspores ~ para_time_point + (1 | clone), data=LA_paraevol) 

overdisp_fun(paraevollogspores)

##       chisq       ratio         rdf           p 
##  44.3092364   0.3142499 141.0000000   1.0000000

summary(paraevollogspores)

## Linear mixed model fit by REML ['lmerMod']
## Formula: logspores ~ para_time_point + (1 | clone)
##    Data: LA_paraevol
## 
## REML criterion at convergence: 259.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.1494 -0.3281  0.1027  0.6344  1.7752 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  clone    (Intercept) 0.02438  0.1562  
##  Residual             0.32184  0.5673  
## Number of obs: 145, groups:  clone, 14
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)       5.02231    0.08202  61.233
## para_time_point3 -0.31407    0.09616  -3.266
## 
## Correlation of Fixed Effects:
##             (Intr)
## par_tm_pnt3 -0.613

drop1(paraevollogspores,test="Chisq")

## Single term deletions
## 
## Model:
## logspores ~ para_time_point + (1 | clone)
##                 npar    AIC    LRT Pr(Chi)   
## <none>               260.88                  
## para_time_point    1 269.32 10.444 0.00123 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#R3 asked about the analysis if we don't log transform the data
paraevolspores <- lmer(spores_tot ~ para_time_point + (1 | clone), data=LA_paraevol) 

overdisp_fun(paraevolspores)

##        chisq        ratio          rdf            p 
## 2.851456e+12 2.022309e+10 1.410000e+02 0.000000e+00

summary(paraevolspores)

## Linear mixed model fit by REML ['lmerMod']
## Formula: spores_tot ~ para_time_point + (1 | clone)
##    Data: LA_paraevol
## 
## REML criterion at convergence: 3813.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.1691 -0.5994 -0.3277  0.3276  4.7299 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  clone    (Intercept) 7.469e+08  27330  
##  Residual             2.042e+10 142903  
## Number of obs: 145, groups:  clone, 14
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)        168349      18911   8.902
## para_time_point3   -54543      24030  -2.270
## 
## Correlation of Fixed Effects:
##             (Intr)
## par_tm_pnt3 -0.658

drop1(paraevolspores,test="Chisq")

## Single term deletions
## 
## Model:
## spores_tot ~ para_time_point + (1 | clone)
##                 npar    AIC   LRT Pr(Chi)  
## <none>               3864.1                
## para_time_point    1 3867.3 5.157 0.02315 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Parasite growth rate

paragrowthparaevol <- ggplot(LA_paraevol_summary, 
                     aes(x=para_time_point, y=meanparagrowth, fill=para_time_point)) + 
  geom_violin(show.legend = FALSE) + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0, seed = 1), alpha = 0.5, show.legend = FALSE) +
  scale_fill_manual(values=c("#C0C0C0", "#f3c483")) + 
  scale_color_manual(values=c("#C0C0C0", "#f3c483")) +
  xlab("Parasite time point") +ylab("Parasite growth rate \n(spores per day)") +
  labs(color= "Parasite time point") + labs(fill="Parasite time point") +
  theme_cowplot()


paragrowthparaevol

paraevolparasitegrowth <- lmer(paragrowth ~ para_time_point + (1 | clone), data=LA_paraevol) 

summary(paraevolparasitegrowth)

## Linear mixed model fit by REML ['lmerMod']
## Formula: paragrowth ~ para_time_point + (1 | clone)
##    Data: LA_paraevol
## 
## REML criterion at convergence: 2778.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3042 -0.6355 -0.2704  0.3463  4.5048 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  clone    (Intercept)   697521  835.2  
##  Residual             14604033 3821.5  
## Number of obs: 145, groups:  clone, 14
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)        5003.8      519.8   9.626
## para_time_point3  -1942.4      644.3  -3.015
## 
## Correlation of Fixed Effects:
##             (Intr)
## par_tm_pnt3 -0.644

drop1(paraevolparasitegrowth,test="Chisq")

## Single term deletions
## 
## Model:
## paragrowth ~ para_time_point + (1 | clone)
##                 npar    AIC    LRT  Pr(Chi)   
## <none>               2815.0                   
## para_time_point    1 2821.9 8.9619 0.002757 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Combined figure

eightpanelplot <- plot_grid(contemplifespanplot, lifespanparaevol, contempreproplot, reproparaevol, contempsporesplot, sporesparaevol, contempparagrowthplot, paragrowthparaevol, labels = "auto", ncol = 2, align = "v")
eightpanelplot

# now add the title
title1 <- ggplot() +                      # Draw ggplot2 plot with text only
  annotate("text",
           x = 0,
           y = 1,
           size = 6,
           label = "Contemporary hosts",
#           hjust = 0.64,
            hjust = 0.5,
           fontface = "bold") + 
  theme_void()

title2 <- ggplot() +                      # Draw ggplot2 plot with text only
  annotate("text",
           x = 0,
           y = 1,
           size = 6,
           label = "Time 3 hosts",
           hjust = 0.5,
           fontface = "bold") + 
  theme_void()

titles <- plot_grid(title1, title2, nrow = 1,   hjust = 0)

titles

eightpanelplot_withtitles <- plot_grid(
  titles, eightpanelplot,
  ncol = 1,
  # rel_heights values control vertical title margins
  rel_heights = c(0.05, 1)
)

eightpanelplot_withtitles

ggsave(here("figures", "eightpanelplot_withtitles.jpg"), eightpanelplot_withtitles, units = "in", width = 8, height = 12, dpi = 300)

Exposed but uninfected animals

Reviewer 3 was asking about exposed but uninfected animals. Let’s make a new dataframe to look at that.

# get rid of time 4 hosts from the full dataset, since we didn't have parasites from that time point
LA_full <- LA_full %>%
  subset(host_time_point != '4')

# filter out ones that are missing spore counts
LA_full <- LA_full %>%
  drop_na(spores_tot)

# Looking for clones that had very few animals in the experiment
LA_full_samplesizes <- LA_full %>%
  count(clone, para_time_point)

# save
write.csv(LA_full_samplesizes, here("tables", "LA_full_samplesizes.csv"), row.names = F)


# Using spore yield to determine if it was infected; this is more sensitive than the 'infection' column, which was based on a quick visual assessment at death -- that inspection missed some that had very low spore densities
LA_full$infected <- ifelse(LA_full$spores_tot == 0, 0, 1)

LA_full$exposureclass <- ifelse(LA_full$para_time_point == "none", 'control',
                                ifelse(LA_full$spores_tot > 0, 'infected', 'exposed'))

LA_full_summary <- LA_full %>%
  group_by(exposureclass, host_time_point, clone, para_time_point) %>%
  summarise(meanlifespan = mean(lifespan), meanclutches = mean(clutches))

## `summarise()` has grouped output by 'exposureclass', 'host_time_point', 'clone'. You can override using the `.groups` argument.

lifespan_exposureclass <- ggplot(LA_full_summary, 
                               aes(x=exposureclass, y=meanlifespan)) + 
  geom_violin() + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, color = '#5c1a33', size = 2) + 
  ylab("Mean lifespan \n(days)")  +
  theme_cowplot() +
  theme(axis.title.x = element_blank())


lifespan_exposureclass

# making a plot comparing exposed vs. unexposed
repro_exposureclass <- ggplot(LA_full_summary, 
                               aes(x=exposureclass, y=meanclutches)) + 
  geom_violin() + 
  geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, color = '#5c1a33', size = 2) +
  ylab("Mean number of clutches")  +
  theme_cowplot() +
  theme(axis.title.x = element_blank())

repro_exposureclass